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FOREWORD 


The  investigation  described  herein  was  performed  by  Dr.  R.  E.  Bunney, 
Division  of  Marine  Resources,  University  of  Washington,  under  Con¬ 
tract  No.  DAAG  17-73-C-0028  for  the  U.S.  A ray  Cold  Regions  Research 
and  Engineering  Laboratory,  sponsored  by  the  Advanced  Research  Project 
Agency  under  ARPA  Order  2096. 

This  contract  was  technically  monitored  by  Dr.  Y.  Nakano,  U.S.  Army 
Cold  Regions  Research  and  Engineering  Laboratory  under  the  instruction 
of  Commander  J.  R.  Seesholtz,  Program  Manager,  ARPA. 
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SUMMARY 


The  principal  object  of  this  experiment  was  to  evaluate  the  feasibil¬ 
ity  of  using  acoustic  pulse  echo  techniques  to  determine  the  subsurface 
structure  of  permafrost.  Laboratory  tests  performed  on  representative 
samples  of  permafrost,  i.e.,  Ottawa  sand,  Hanover  silt  and  Goodrich  clay, 
show  that  such  acoustic  methods  are  reasonable,  provided  the  acoustic  wave¬ 
lengths  are  long  compared  to  the  radius  of  the  scattering  centers  and  to 
the  transition  regions  between  adjacent  constituents  (e.g.,  sand,  silt, 
clay,  ice)  of  the  permafrost. 

Our  velocity  measurements  showed  that  either  previous  measurements 
have  been  significantly  high  or  the  saturation  levels  of  our  samples  were 
lower  than  anticipated.  The  test  sections  deteriorated  visibly  during  the 
experimentation  process.  There  was  obviously  dehydration  due  to  sublima¬ 
tion  of  the  surface  moisture,  but  the  extent  of  dehydration  in  the  in¬ 
terior  of  the  sample  is  not  known.  The  velocity  measurements  suggest 
saturation  levels  of  35-40%. 

Acoustic  attenuation  measurements  were  performed  and  the  results  are 
compared  with  scattering  theory.  The  additional  attenuation  expected  from 
dissipative  processes  is  discussed  but  could  not  be  evaluated  because  the 
dimensions  of  our  samples  were  not  large  enough  to  allow  the  long  wave¬ 
length  tests  necessary  to  determine  the  parameters  required  in  the  theory. 
The  results  of  the  tests  did  indicate,  however,  that  if  significant  depth 
is  to  be  obtained  frequencies  of  less  than  10  kHz  must  be  utilized.  It 
is  also  clear  that  an  efficient  method  of  coupling  the  energy  into  the  me¬ 
dium  must  be  investigated.  From  the  results  of  our  theories  and  measure¬ 
ments,  there  appears  to  be  no  significant  difference  between  using  the 
compressional  and  using  the  transverse  acoustic  mode. 

The  reflectivity  at  permafrost-ice  interfaces  was  investigated  and 
the  results  are  compared  with  theory.  Although  both  the  predicted  and 
measured  losses  in  the  reflected  wave  are  appreciable,  they  do  not  appear 
to  preclude  using  acoustic  pulse  echo  techniques,  provided  the  wavelengths 
are  long  enough  to  reduce  the  attenuation  to  acceptable  levels.  Extension 
of  the  theory  to  include  transition  layers  between  constituents  of  the 
medium  is  also  discussed. 


INTRODUCTION 


Many  operations  in  arctic  regions  require  detailed  analysis  of  the 
shallow  (<50  ft)  permafrost  subsurface.  This  information  is  particularly 
important,  for  example,  in  the  construction  of  oil  pipelines;  in  such 
cases,  however,  the  survey  of  the  permafrost  subsurface  should  not  only 
be  efficient  and  not  too  costly  but  should  have  a  minimal  environmental 
impact.  Determining  the  feasibility  of  using  acoustic  pulse  echo  methods 
to  locate,  and  possibly  evaluate,  boundaries  and  discontinuities  in  the 
potmafrost  has  been  the  subject  of  the  study  reported  here. 

The  ability  to  use  acoustic  ranging  to  locate  inclusions  and  dis¬ 
continuities  in  permafrost  formed  by  dry  frost  regions,  the  subpermafrost 
boundary,  refrozen  regions  and/or  upcroppings  of  base  rock  (all  of  which 
may  serve  as  foundations  for  substantial  structures)  or  to  locate  ice 
lenses,  taliks,  etc.  (which  can  not)  is  dependent  upon  the  acoustic  prop¬ 
erties  of  the  respective  media  and  the  transition  zones  between  them. 
Specifically,  the  densities  and  either  the  longitudinal  and  transverse 
acoustic  velocities  or  the  elastic  moduli  of  the  permafrost  constituents 
must  be  very  well  known  to  determine  the  reflectivity  at  the  interfaces; 
the  attenuation  as  a  function  of  frequency  must  also  be  known  to  determine 
the  depth  to  which  the  method  will  be  reliable. 

A  great  deal  of  work  has  been  done  on  measuring  the  velocity  of 
sound  in  permafrost  using  seismic  techniques  (see  the  review  by  Barnes1 
and  Zykov*) .  Hunter3  has  recently  applied  seismic  techniques  to  mapping 
subsurface  structures  and  has  also  reported  the  results  of  in  situ  com- 
pressional  and  transverse  wave  velocity  measurements.  AlthougH- labora¬ 
tory  measurements  of  the  acoustic  velocity  in  various  constituents  of 
permafrost  have  been  performed  for  some  time1*-*  and  many  theories7”13 
that  are  applicable  to  the  propagation  of  sound  in  permafrost  have  been 
set  forth,  only  in  the  last  few  years  has  significant  progress  been  re¬ 
ported  on  determining  the  sound  propagation  properties  of  frozen  soils. 
Kapler1*  used  a  resonant  method  to  determine  the  elastic  moduli  of  sev¬ 
eral  representative  samples  of  permafrost  and  ice  as  a  function  of  tem¬ 
perature.  Timur13  discusses  the  compressional  velocity  in  terms  of  the/ 
ice-water-air  constituents  of  the  media  and  relates  the  observed  increase 
in  velocity  to  the  cementation  of  the  sand  grains  with  ice.  He  also^ob- 
served  a  hysteresis  effect  in  the  velocity  coinciding  with  the  freeze- 
thaw  cycle  and  related  it  to  interfacial  forces  and  a  slight  salinity  of 
the  interstitial  liquids.  Nakano  and  his  co-workers13 ”13  have  performed 
a  series  of  careful  laboratory  experiments  on  samples  of  Ottawa  sand, 
Goodrich  clay  and  Hanover  silt,  measuring  the  compressional.  velocity  using 
pulse  echo  techniques  and  measuring  the  shear  velocity  using  Snell's  law 
of  critical  angle  methods.  For  the  case  of  Ottawa  sand,  the  velocity  was 
related  not  only  to  temperature,  where  the  hysteresis  effect  of  Timur  was 
observed,  but  also  to  the  saturation  of  the  sample. 


As  discussed  previously,  the  attenuation  of  the  sound  as  a  function 
of  frequency  must  be  well  known  if  the  investigator  is  to  maximize  both 
the  depth  of  penetration  into  the  medium  and  the  resolution  with  which 
the  subsurface  inclusions  and  discontinuities  are  observed.  Unlike  ve¬ 
locity,  there  unfortunately  appears  to  have  been  little  work  performed 
in  this  area,  at  least  in  the  laboratory.  Anomalously  high  attenuation 
has  been  observed  by  Hunter3  and  by  Press  and  Dobrin19  during  seismic  ex¬ 
periments  in  regions  of  thin  permafrost.  However,  since  Donato20  has 
shown  that  for  this  condition  the  attenuation  is  related  to  the  ratio  of 
the  wavelength  to  the  thickness  of  the  sample,  this  anomaly  is  probably 
related  to  Lamb-like  propagation  modes21  in  the  medium  rather  than  its 
gross  properties.  An  excellent  review  article  (with  extensive  bibliog¬ 
raphy)  on  the  attenuation  of  acoustic  waves  at  seismic  frequencies 
has  been  written  by  Jackson  and  Anderson.22  However,  because  of  the  long 
wavelengths  used  in  seismic  work,  the  results  in  those  studies  may  not  be 
applicable  to  the  problem  being  studied  here  because  of  the  resolution 
that  is  required  to  adequately  analyze  the  subsurface  structure. 
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EXPERIMENTS 

Samples 

The  data  reported  were  obtained  on  idealized,  i.e.,  homogeneous  and 
isotropic,  samples  of  representative  permafrost  constituents  furnished 
by  the  U.S.  Corps  of  Engineers'  Cold  Regions  Research  and  Engineering 
Laboratory.  Specifically,  we  used  samples  of  20-30  Ottawa  sand,  Hanover 
silt  and  Goodrich  clay,  formed  in  cylinders  approximately  6  inches  in 
diameter  and  18  inches  long.  To  prepare  the  samples,  the  dry  soil  was 
saturated  with  distilled  water,  and  then  placed  in  a  mold.  The  mold  was 
placed  in  a  cold  room  and  subjected  to  rotary  motion  to  prevent  localized 
concentrations  due  to  gravity  while  the  sample  was  freezing.17 

During  the  course  of  the  experiment,  some  evidence  of  deterioration 
became  apparent.  Cross  sections  of  the  sample  indicate  this  deteriora¬ 
tion  certainly  was  caused  by  sublimation  at  the  surface  and  perhaps  con¬ 
tributed  to  a  general  degeneration  of  the  sample.  It  most  assuredly  was 
a  contributing  factor  to  the  anomalous  results  that  are  discussed  later. 

Results  of  Velocity  Measurements 

Figure  1  is  a  block  diagram  of  the  instrumentation  used  to  obtain 
the  velocity,  attenuation  and  reflectivity  measurements.  In  this  system 
a  pulse  timing  generator  furnishes  pulse  length  and  rate  information  to 
a  pulsed  power  oscillator,  and  a  trigger  signal  to  the  oscilloscope.  The 
oscillator,  in  turn,  provides  a  cw  pulse  of  the  desired  frequency,  rate  and 
length  to  an  acoustic  transducer  coupled  to  the  sample.  The  transmitted 
acoustic  signal  is  received  by  a  second  transducer  fixed  to  the  opposite 
side  of  the  test  section,  amplified,  and  transmitted  to  the  oscilloscope. 
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Figure  1.  Block  diagram  of  instrumentation. 


The  time  of  flight  of  the  acoustic  pulse  in  the  sample  is  determined 
by  measuring,  with  a  calibrated  time  delay,  the  oscilloscope  sweep  time 
from  the  origin  to  the  first  received  signal  and  to  the  subsequent  arrivals 
that  represent  multiple  reflections  in  the  test  section.  Typical  results 
are  shown  in  Figures  2  and  3. 

On  relatively  homogeneous  materials,  measurements  using  this  tech¬ 
nique  are  normally  accurate  to  within  1  1/2%,  with  the  largest  error  in 
the  measurement  of  length.  For  the  materials  being  discussed  here,  how¬ 
ever,  the  velocity  measurements  are  not  that  accurate  because  of  loss  of 
resolution  due  to  multipath  spreading  of  the  pulse  during  transmissions 
and,  in  some  cases,  because  of  the  long  wavelengths  required  due  to  the 
high  attenuation. 

The  results  of  the  velocity  measurements  are  shown  in  Figures  4-6 
and  the  mean  values  are  tabulated  in  Table  I.  During  the  tests,  the 
sample  temperature  was  maintained  as  nearly  constant  as  possible  at  -15°C. 


Table  I.  Tabulation  of  average  velocities . 


Material _ Ottawa  Sand _ Hanover  Silt _ Goodrich  Clay 


Comp.  Vel.  3417  1  435  m/sec  3096  1  309  m/sec 

Trans.  Vel.  2266  ±  424  m/sec  1805  ±  335  m/sec  2065  ±  206  m/sec 
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Figure  2.  Typical  result  of  velocity  experiment  on  Ottawa  sand . 
Frequency  =  64  kHz,  sweep  =  SO  psec/cm. 


Figure  4.  Velocity  va  frequency ,  Ottawa  eand. 
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Figure  S.  Velocity  os  frequency,  Banover  eilt 
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The  results  differ  appreciably  from  those  of  Nakano  and  Froula1*  (Figures 
7  and  8)  at  the  same  temperature.  Although  the  cause  of  the  disparities 
is  not  known,  dehydration  of  the  samples  is  a  strong  possibility.  This 
hypothesis  is  substantiated  by  Nakano's  and  Arnold's17  measurements  of 
the  velocity  in  frozen  Ottawa  sand  vs  saturation  and  temperature  (see 
Figures  9  and  10,  and  Table  II).  These  data  indicate  that  our  Ottawa 
sand  samples,  rather  than  being  totally  saturated,  actually  were  only 
saturated  to  35-40%.  Unfortunately,  similar  data  are  unavailable  for 
Hanover  silt  and  Goodrich  clay.  These  materials  would  be  expected,  how¬ 
ever,  to  behave  in  a  fashion  similar  to  that  of  Ottawa  sand.  Thus,  the 
lower  velocity  values  observed,  particularly  in  the  case  of  Hanover  silt, 
are  probably  due  to  a  lower  saturation  level  and  are  not  necessarily 
surprising. 


Table  II.  Description  of  Specimens 


Specimen 

Thickness  d 
(cm) 

Net 

Density 

(g/cm1) 

Dry 

Density 

(g/cm') 

Nater 

Content* 

C%) 

Voi Patio 

Ice 

Saturationf 

(%) 

1 

2.394 

1.962 

1.609 

21.91 

0.6470 

99.70 

2 

2.373 

1.863 

1.598 

16.59 

0.6583 

74.19 

3 

2.437 

1.773 

1.554 

14.09 

0.7053 

58.82 

4 

2.399 

1.771 

1.589 

11.47 

0.6677 

50.57 

S 

2.363 

1.71S 

1.56S 

9.59 

0.6933 

40.72 

6 

2.400 

1.717 

1.S92 

7.85 

0.6646 

34.77 

7 

1.321 

1.670 

1.590 

5.04 

0.6667 

22.27 

8 

1.339 

1.604 

1.562 

2.66 

0.6965 

11.24 

9 

1.318 

1 .588 

1.572 

1.02 

0.68S8 

4.40 

10 

• 

1.311 

1.580 

1.580 

0.00 

0.6772 

0.00 

Height  percent  of  water  referred  to  the  unit  weight  of  dry  soil. 

^Weight  percent  of  ice  referred  to  the  weight  of  ice,  which  could  fill  total  void, 
(after  Nakano  and  Arnold17) 


Results  of  Attenuation  Measurements 


The  attenuation  of  the  acoustic  wave  in  the  medium  is  a  principal 
concern  in  determining  both  the  resolution  and  depth  that  can  be  attained 
with  pulse  echo  techniques.  This  attenuation,  whether  originating  from 
dissipative  considerations  (i.e.,  internal  friction,  thermal  dissipation, 
viscous  slippage  at  crystal  boundaries,  etc.)  or  from  acoustic  scattering 
from  inhomogeneities  in  the  medium,  is  very  sensitive  to  frequency.  In 
fact,  the  resolution  of  the  technique  is  inversely  proportional  to  the 
frequency  of  the  acoustic  wave.  Thus  to  optimize  his  experiments,  the 
investigator  must  be  able  to  predict  the  attenuation  as  a  function  of 
frequency. 
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VELOCITY  ( 


Hgur e  The  dilatational  wave  velocity  va  temperature  for:  A)  20-30 

Ottawa  Band ,  wet  density  p  »  2.20  g/am* ,  B)  Hanover  siltt 
p  »  2.03  g/am*t  and  C)  Goodrich  clay,  p  *  1.80  g/arn *  under 
fully  wdter-eaturated  conditions,  notice  the  discrepancy 
between  the  freezing  and  thawing  cycles  for  B  end  C.  Ho 
discrepancy  was  observed  for  Ottawa  sand,  (after  Nakano1*) 


Figure  9. 


Figure  10. 


TOPVMTUNi  PC) 


Dilatational  wave  velocity  versus  temperature 
under  various  conditions  of  ice  saturation. 
The  number  assigned  to  each  curve  corresponds 
to  the  specimen  number,  (after  Nakano  and 
Arnold17) 


Shear  wave  velocity  versus  temperature  under 
various  conditions  of  ice  saturation.  The 
umber  assigned  to  each  curve  corresponds  to 
the  specimen  number,  (after  Nakano  and 
Arnold17) 
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An  acoustic  wave  propagating  in  a  nonhomogeneous  medium  will  inter* 
act  with  the  inhomogeneities  and  part  of  the  energy  will  be  scattered  out 
of  the  sound  beam.  When  the  wavelength  of  the  sound  is  large  compared  to 
the  size  of  these  scattering  centers,  the  interaction  is  minimal  and  the 
energy  loss  is  small  compared  to  that  caused  by  classical  dissipation. 

For  this  situation,  the  acoustic  attenuation  can  be  predicted  by  the 
long  wavelength  theories  discussed  in  Appendix  A.  However,  when  the 
size  of  the  scattering  center  is  a  significant  portion  of  the  wavelength, 
the  attenuation  due  to  scattering  cannot  be  neglected. 

To  determine  the  attenuation  due  to  scattering,  consider  an  isotropic 
solid  elastic  medium  containing  a  concentration  of  scatterers.  The  inhomo- 
geneities  in  the  medium  can  be  mathematically  characterized  by  a  generalized 
function. 


D(S,I,a,r. . .)  -  SCsiJICsi.s^GCa,?)  ...,  (1) 

where  S(si)  is  a  structure  distribution  function  which  describes  the 
relative  importance  of  the  various  types  of  scatterers  in  the  medium, 
I(si.sj)  is  a  distribution  function  which  describes  the  interaction 
between  scatterers  and  multiple  scattering  of  the  sound,  and  G(a,r)  de¬ 
scribes  the  distribution  of  the  size  and  position,  etc.  of  the  scatterers. 

Assuming  a  plane  wave  moving  in  the  positive  x-direction,  the  intensity 
of  .the  sound  is  given  by 


I  -  ID  e"20*  ,  (2) 

where  a  is  the  attenuation  constant,  x  is  the  penetration  distance  into 
the  material  and  IQ  is  the  intensity  of  the  original  sound  wave.  If  the 
energy  losses  are  all  considered  to  be  caused  by  scattering  processes, 

Eq.  (2)  may  be  rewritten 

i  -  i0  ®'r*  .  (S) 

where  T  is  the  generalized  scattering  cross  section  and  is  governed  by  the 
size,  shape,  position,  etc.  of  the  scatterers  in  the  medium.  In  terms  of 
Eq.  (1),  r  may  be  written 


T  - 


//...D(S,I,G.r,...)  y,  (a)drda... 
r  a 


(4) 


where  now  y^(a)  is  the  scattering  cross  section  of  a  single  independent 
scatterer.  The  attenuation  is  then  given  by 
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This  equation  states  that  the  attenuation  is  determined  by  the  scat¬ 
tering  cross  section  and  the  distribution  of  the  scatterers.  The  cross 
section  is  dependent  upon  the  interaction  of  the  wave  with  each  individual 
scatterer,  while  the  distribution  function  is  the  probability  of  the  wave 
encountering  that  scatterer  in  the  medium. 

For  purposes  of  calculation,  permafrost  is  considered  to  be  composed 
of  an  array  of  elastic  scatterers,  i.e.,  granules  of  Ottawa  sand,  etc., 
imbedded  in  an  ice  matrix.  To  calculate  the  scattering  cross  section  of 
a  single  scatterer,  the  problem  shown  in  Figure  11  must  be  resolved. 
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WAVE 
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Figure  11.  Representation  of  the  problem. 


Incorporating  the  methods  of  Ying  and  Truell23  and  of  Corbato"  and 
Uretsky,24  and  using  the  experimentally  measured  values  of  velocity  and 
density  as  reported  by  Nakano,18  the  product  of  the  scattering  cross 
section  and  the  square  of  the  frequency  (yv2)  were  calculated  numerically 
as  a  function  of  the  circumference-to-wavelength  ratio  for  both  a  shear 
and  compressional  incident  plane  wave.  The  results  of  this  calculation 
are  shown  in  Figures  12-14.  In  the  frequency  region  for  ka«l,  i.e., 
for  wavelengths  much  larger  than  the  scatterer  radius,  the  scattering 
cross  section  is  of  the  form 
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-  v' 


(6) 


[\ 


which  is  the  result  expected  from  standard  Rayleigh28  scattering.  In 
general,  the  scattering  cross  sections  for  ka<l  and  ka>15  can  be 
approximated  by 
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Hanover  silt 
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r shear 


7.86  x  10“9  v2*41  a2*41  +  6.87  x  10" 
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dB/cm2 


where  the  frequency  is  in  megahertz  and  the  scatterer  radius,  a,  is  in 
centimeters . 

To  evaluate  the  attenuation  due  to  scattering  from  these  expressions, 
the  distribution  functions  in  the  medium  must  be  either  known  or  assumed. 
Because  the  necessary  data  are  presently  unavailable,  and  to  minimize 
the  complexity  of  the  calculation,  the  following  assumptions  have  been 
made. 

1.  No  interaction  exists  between  the  scattering  centers,  and  the 
wave,  once  scattered,  will  not  be  re-scattered  back  into  the 
sound  beam.  This  approximation  is  not  extremely  valid  for 
materials  with  as  high  a  concentration  of  scattering  centers 
as  those  being  discussed  here  and  will  tend  to  make  the  pre¬ 
dictions  somewhat  higher  than  if  multiple  scattering  had  been 
considered. 

2.  There  are  N  scattering  centers  per  unit  volume  and  the  centers 
are  uniformly  distributed  throughout  the  medium.  With  the 
exception  of  the  Goodrich  clay  sample,  this  appears  to  be  a 
reasonable  assumption.  It  will,  of  course,  not  be  valid  for 
in  situ  measurements . 

3.  All  of  the  scattering  centers  are  of  the  same  size,  with  a 
radius  a0.  Since  the  materials  used  in  the  manufacture  of  the 
test  samples  were  filtered  for  size,  this  assumption  is  per¬ 
fectly  reasonable.  During  field  tests,  this  assumption  becomes 
invalid.  However,  if  radius  a<j  is  selected  as  the  average 
scatterer  size,  then  it  can  be  argued  that  only  the  larger 
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inhoinogeneities  in  the  distribution  will  affect  the  scatter¬ 
ing.  This  assumption  will  make  the  calculated  attenuation 
less  than  that  expected  during  the  experiment. 


Using,  these  assumptions,  the  attenuation  due  to  scattering  may  be 
written 


Scattering  *  2  ^ao^  »  (10) 

which  predicts  the  "first  order"  attenuation  to  be  expected  from  the 
scattering  mechanism.  Using  the  curves  of  Figures  12  through  14  and 
the  relation 


ka 


2ttv 

C  a  » 


(ID 


where  v  is  the  frequency  in  hertz,  C  is  the  velocity  in  centimeters  per 
second,  and  a  is  the  radius  of  the  scatterers  in  centimeters,  the  atten¬ 
uation  due  to  scattering  for  values  ka<l  and  ka>15  can  be  approximated 
by: 
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ccmp  5.50  x  10“2  v2-*  ♦  1.15  x  10-* 

dB/m 


dB/m 


shear 


1.14  x  10'2  va  ♦  1.03  x  10 


(12) 


(13) 
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where  v  is  in  megahertz.  For  values  of  l<ka<  IS,  the  approximations  of 
Eqs.  (12)  through  (14)  are  still  more  or  less  valid  for  order  of  magnitude 
calculations,  but  will  invariably  yield  predictions  that  are  significantly 
low. 


Equations  (12)  through  (14)  have  been  solved  as  a  function  of  frequency 
and  the  solutions  are  shown  in  Figures  15  through  17.  For  the  frequencies 
considered  during  these  experiments,  relatively  good  agreement  was  ob¬ 
served  between  the  theory  and  the  experiment,  with  the  exception  of  the 
Goodrich  clay  samples.  The  reason  for  the  variance  between  the  theory 
and  the  results  of  the  Goodrich  clay  experiment  can  be  seen  in  Figure  18. 
The  theories  assume  both  total  saturation  and  homogeneity.  Earlier  it 
was  speculated  that  during  the  course  of  the  experiment  the  test  samples 
degenerated  appreciably  because  of  dehydration.  This  degeneration  was 
assumed  to  be  limited  to  sublimation  of  the  surface  areas,  but  subsequent 
analysis  of  the  velocity  data  has  shown  that  the  samples  were  not  totally 
saturated  at  the  time  of  the  experiment.  Figure  18  indicates  that  homo¬ 
geneity  also  was  not  necessarily  achieved  in  the  test  samples. 

At  the  beginning  of  this  section  we  mentioned  the  attenuation  loss 
due  to  dissipative  mechanisms,  i.e.,  thermal  losses,  internal  friction, 
viscous  slippage,  etc.  These  processes  have  not  been  specifically  eval¬ 
uated  for  the  permafrost  specimens  studied  here,  primarily  because  the 
dimensions  of  the  test  samples  precluded  using  acoustic  frequencies  that 
were  very  long  compared  to  the  radius  of  the  scattering  centers.  Further 
experimentation  should  be  performed  to  acquire  the  long-wave length  param¬ 
eters  necessary  to  evaluate  the  long  wavelength  solutions  discussed  in 
Appendix  A. 

Results  of  Reflectivity  Measurements 

The  energy  of  an  acoustic  wave  encountering  a  boundary  between  two 
elastic  media  is  divided  between  reflected  and  transmitted  modes  (see 
Figure  19).  To  determine  the  feasibility  of  analyzing  the  subsurface 
structure  of  permafrost,  the  relative  amount  of  energy  in  the  reflected 
■odes  must  be  known.  In  Appendix  A,  the  problem  of  two  ideal  elastic 
media  in  rigid  contact,  i.e.,  presuming  no  slippage  at  the  boundary,  has 
been  solved  as  a  function  of  incident  angle  for  both  compressional  and 
transverse  acoustic  waves  incident  on  a  permafrost- to- ice  boundary.  The 
results  of  this  calculation  are  shown  in  Figures  20-25  and  indicate  that, 
for  normal  incidence,  a  minimum  loss  of  7  to  11  dB  is  to  be  expected  at 
the  boundary,  depending  on  the  mode  of  the  incident  wave  and  the  perma- 
frost  material;  for  an  Ottawa  sand-Goodrich  clay  interface  (Figures  26 
and  27),  losses  exceeding  13  dB  can  be  anticipated  in  the  reflected  wave. 
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ATTENUATION  (dB/m) 


FREQUENCY  (Hr) 


Figure  16.  Attenuation  ve  frequency  (eoattering  only),  Hanover  eilt. 


22 


Figure  18.  Photograph  of  aro88  section  of  Goodrich  clay  sample. 


Relative  energy  distribution  of  aoouatio  nave  incident  on  permafroet~ioe 
Ottawa  eand-ioe  (oompreeeional  wave  incident). 


iva  energy  distribution  of  aeoustie  wave  incident  on 
ter  eilt-ioe  (oampreeeiontxl  wave  incident). 


t tion  of  aoouetio  wave  incident  on 


Relative  energy  distribution  of  aoouetio  wave  incident  on 
Goodrich  clay-ice  (oompreeeional  wane  incident). 


Figure  25.  Relative  energy  distribution  of  acoustic  wave  incident  on  permf'r oat-ice 
Goodrich  clay-ice  (traneveree  wave  incident). 


Figure  27.  Relative  energy  distribution  of  aoouetio  wave  incident  on  boundary  of  two  permafrost  layers, 
Ottawa  sand-Goodrioh  olau  (transverse  wave  incident). 
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To  experimentally  measure  the  reflectivity  theoretically  predicted 
above,  the  experiment  diagrammed  in  Figure  28  was  performed.  In  this 
test,  an  ice-permafrost  boundary  was  formed  as  shown  in  Figure  29.  A 
hydrophone  was  imbedded  in  the  ice  a  known  distance  from  both  the  surface 
and  the  interface.  An  acoustic  transducer  was  coupled  to  the  ice  surface 
and  pulsed.  The  amplitude  of  the  acoustic  wave  was  then  recorded  by  the 


hydrophone  both  for  the  initial  pulse  and  for  the  reflection  from  the 
interface.  A  typical  result  is  shown  in  Figure  30.  Correcting  for  at¬ 
tenuation  losses  in  the  ice,  it  was  then  possible  to  determine  the  energy 
loss  at  the  boundary.  The  results  of  this  measurement  are  shown  in  Table 
III.  With  the  exception  of  Hanover  silt,  the  losses  observed  during  the 
test  are  clearly  significantly  larger  than  the  theory  predicted.  There 
are  many  reasons  for  this  result,  including  irregularities  at  the  inter¬ 
face,  inhomogeneous  media,  etc.,  but  not  the  least  of  them  is  the  exist¬ 
ence  of  a  transition  layer.  Figure  31  is  a  diagram  of  the  transition  layer 
problem. 

Table  III .  Results  of  reflectivity  measurements 

(energy  loss  at  the  permafrost-ioe  boundary). 


Material 

Predicted 

Value 

Measured 

Value 

Ottawa  Sand  -  Ice 

7.0  dB 

12.9  dB 

Hanover  Silt  -  Ice 

9.9  dB 

10. 5  dB 

Goodrich  Clay  -  Ice 

8.3  dB 

16.2  dB 
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Figure  29.  Photograph  of  test  sample  used  in  reflectivity  measure¬ 
ment  for  ice-clay  interface.  (Note  the  condition  of 
the  clay  surface.) 


35 


Figure  30.  Typical  result  of  reflectivity  measurement. 

(Arrow  1  denotes  initial  pulse;  arrow  2  de¬ 
notes  reflected  pulse.) 


Figure  31.  Diagram  of  transition  layer  problem. 


The  reflectivity  of  an  acoustic  plane  wave  normally  incident  on  a 
layer  between  two  isotropic  media  for  which  the  velocity  and  density  are 
dependent  upon  the  depth  into  the  layer  is  theoretically  considered  in 
Appendix  A  for  the  case  of  C  =  (l+ax)n  as  a  function  of  the  thickness- 
to-wavelength  ratio  in  the  incident  medium,  k0h.  The  results  of  this  cal¬ 
culation  for  the  special  cases  of  n=l  (a  linear  transition  between  media) 
are  shown  in  Figures  32-34.  This  calculation  demonstrates  that  the  ef¬ 
fect  of  the  layer  is  insignificant  for  wavelengths  that  are  long  compared 
to  the  thickness.  However,  when  the  wavelength  becomes  on  the  order  of 
the  thickness,  the  reflectivity  becomes  oscillatory  with  increasing  kjjh. 
Thus,  some  knowledge  of  the  interface  between  permafrost  materials  will 
be  required  to  allow  the  experimenter  to  use  acoustic  frequencies  that 
will  result  in  the  maximum  energy  possible  being  reflected  from  the 
boundary . 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  feasibility  of  utilizing  acoustic  range-echo  techniques  to  effi¬ 
ciently  analyze  subpermafrost  structure  to  significant  depths  depends  upon 
several  factors.  The  more  critical  problems,  those  associated  with  the 
velocity  and  attenuation  of  the  acoustic  wave  in  the  medium,  with  the 
acoustic  impedance  mismatch  at  boundaries  between  various  constituents  of 
permafrost,  and  with  the  modulation  caused  by  possible  transition  layers, 
have  been  the  subject  of  this  study.  The  results  of  the  research  indicate 
that,  if  sufficient  coupling  between  the  sound  source  and  the  medium  can  be 
obtained,  and  if  the  acoustic  wavelength  is  long  compared  to  the  majority 
of  the  scattering  centers  and  to  the  transition  layer  between  media,  the 
utilization  of  pulse  echo  techniques  could  be  a  reasonable  solution  to  the 
analysis  problem. 

Unfortunately,  the  results  of  the  study  are  not  as  conclusive  as  one 
would  like.  It  is  therefore  recommended  that  further  basic  research  on 
the  acoustic  properties  of  permafrost  materials  be  performed.  Specific¬ 
ally,  further  test  data  should  be  obtained  on  the  acoustic  velocity  and 
attenuation  of  permafrost  samples  for  which  the  physical  properties,  i.e., 
saturation,  homogeneity,  etc.,  have  been  carefully  controlled.  Once  the 
acoustic  properties  of  these  "idealized"  permafrost  samples  are  well  under¬ 
stood,  the  tests  should  be  expanded  to  include  studying  the  acoustic  prop¬ 
erties  of  cores  extracted  from  the  environment,  provided  techniques  exist 
to  acquire  these  samples  without  disturbing  their  physical  composition. 

If  not,  in  situ  measurements  should  be  performed  prior  to  the  extraction 
of  the  permafrost  cores. 

Finally,  it  is  of  great  importance  to  initiate  a  program  to  study 
the  physical  properties  of  transition  layers  between  various  constituents 
of  permafrost.  Although  the  experimenter  can  optimize  his  tests  for  depth 
and  resolution  by  adjusting  the  acoustic  frequency,  unless  che  transition 
zones  can  be  successfully  modeled,  successful  and  efficient  subsurface 
permafrost  reflectivity  measurements  will  never  be  achieved. 
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Figure  83.  Reflectivity  from  transition  layer ;  Hanover  silt -ice;  velocity  m  (1+AX) 


0*175  f-  □  - COMPRESSIONS  INCIDENT 


Figure  84.  Reflectivity  from  transition  layer;  Goodrich  olay  -  ice;  velocity  »  ( 1+AI) 


APPENDIX  A 


THEORIES 


I.  THEORY  OF  ACOUSTIC  REFLECTION  AND  TRANSMISSION  AT  AN  INTERFACE 
A.  Ideal  Case 

Consider  an  acoustic  plane  wave  striking  an  interface  between  two 
elastic  media  as  shown  in  Figure  Al.  It  is  assumed  that  each  medium  is 
isotropic  and  homogeneous.  Complexities,  such  as  viscoelastic  properties, 
are  not  considered.  Each  medium  is  characterized  by  its  density.  Lame' 
constant  and  shear  modulus.  Finally,  it  is  assumed  that  strong  coupling 
exists  at  the  interface. 


Figure  Al.  Accrue  tio  wave  reflection  and  tremmieeion  at  an  interface. 


The  displacement  potential  functions  are  then  given  by 


42 


♦j  -  [6L.-l00It.  A2 .‘v| 


where 


and 


♦j  ■  .‘(“-M  [(i-4l)  ^“S*.  Bj.10®1] 

♦2  .  .““-M  A'.-ioL* 

*  .  BviaI*  _ 

11  incident  compressional  wave 

0  incident  transverse  wave 

3  »  kQ  sin  0Q  *  kg  sin  0g  »  kL  sin  0^  «  1^  sin  0T 

ao  "  kocos0o 

°S  *  kScos0S 

°L  *  kLCOS  ®L 

(Xj,  *  ky  cos  0T 


(Al) 


At  the  boundary  x»0,  the  normal  displacement,  pressure  and  shear  stress 
are  continuous.  Also,  since  strong  coupling  has  been  assumed,  the  tangen- 


tial  displacement  is  also  continuous.  Mathematically,  this  is  expressed 
as 
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where 


u  »  grad  <p  *  curl  $ 

I  »  (0,0,i|0 


10  i« 


CA3) 


and  i  and  j  denote  direction.  Applying  the  boundary  conditions  then  gives 
four  simultaneous  equations 


which  have  been  solved  using  program  SOLID  in  Appendix  B  to  calculate  the 
energies  of  reflection  and  transmission  at  the  interface  for  several  ex¬ 
amples.  The  results  of  the  calculation  are  shown  in  Figures  20-27  of  the 
text . 

B.  Extension  of  the  Theory  to  Include  a  Transition  Layer 

The  problem  to  be  solved  is  diagrammed  in  Figure  A2.  An  incident 
plane  compressional*  wave  is  incident  on  the  boundary  of  a  layer  in  which 
the  velocity  and/or  the  density  is  a  function  of  depth  into  the  medium. 
The  displacement  potentials  are  given  by 

♦0  -  C*il‘°x  *  A 2  «ik»x)  tiM 

♦j  -  (Bj  «‘iktx  ♦  B2  eiklX)  e1"*  (AS) 
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“Development  of  the  theory  for  an  incident  transverse  wave  is  equivalent 
to  the  development  presented  here. 


(A4) 
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Figure  A2.  Diagram  of  transition  layer  problem. 

where  now  kL  is  a  function  of  x.  At  the  boundaries  x=0  and  x*h  the  normal 
displacement  and  pressure  are  continuous,*  where  the  displacements  and 
stresses  are  defined  by  Eq.  (A3).  Evaluating  the  boundary  conditions 
yields  four  simultaneous  equations  to  be  solved  for  the  amplitudes: 

v1-^  -  kL|  £BrB2> 

1  XmO 

9k 

p0a-*2)  -fit!  »i-b2> 

l  L  x*o  i 

(vh  ^r)(Bi*'lkLh-B2*lkLl’)  |x.h  ■  kc‘*'lkch  <A6) 


*For  the  case  of  a  shear  wave  incident,  the  tangential  displacement 
and  shear  stress  are  continuous  where  it  has  been  assumed  that  strong 
coupling  exists  at  the  boundaries. 


which  are  completely  general. 

Consider  the  special  case  for  which  the  functional  dependence  of  the 
velocity  in  the  transition  region  is  given  by 

CL  -  CQ  (1+ouc) "  ,  (A7) 

where  a  is  the  gradient  and  n  is  the  order  of  dependence  and  need  not 
necessarily  be  an  integer.  It  can  then  be  shown  that  the  amplitude  of 
the  reflected  wave  is  given  by 
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These  equations  have  been  evaluated  using  the  computer  program  TRNTION 
given  in  Appendix  B  for  the  special  cases  of  n=l  and  n=2.  The  results 
of  this  calculation  are  shown  in  Figures  32-34  of  the  text. 

II.  APPROXIMATION  OF  ATTENUATION 

The  attenuation  of  a  sound  wave  propagating  in  a  material  that  is  not 
perfectly  elastic,  homogeneous  and  isotropic  can  be  considered  to  be  the 
result  of  two  mechanisms:  (1)  dissipation  processes  originating  from  in¬ 
ternal  friction,  anelastic  behavior  of  the  material,  thermal  dissipation, 
viscous  slippage  at  crystal  boundaries,  etc.,  and  (2)  scattering  originat¬ 
ing  from  the  interaction  of  the  acoustic  wave  with  scattering  centers  in 
the  medium.  The  relative  contribution  of  each  of  these  mechanisms  to  the 
total  attenuation  depends  on  the  frequency  of  the  sound  wave.  At  low  fre¬ 
quencies,  the  wavelength  of  the  sound  is  very  large  compared  to  the  scat¬ 
tering  centers;  therefore,  the  scattering  cross  section,  i.e.,  the  rela¬ 
tive  amount  of  energy  scattered  out  of  the  incident  wave,  is  extremely 
small  and  the  attenuation  is  due  almost  entirely  to  dissipation.  As  the 
frequency  increases,  the  attenuation  due  to  scattering  becomes  more  im¬ 
portant  until,  when  the  wavelength  becomes  approximately  on  the  order  of 
the  size  of  the  scattering  center,  scattering  predominates.  Thus,  to 
assess  the  frequency  dependence  of  the  attenuation,  the  dissipation  and 
scattering  processes  must  both  be  evaluated. 
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A.  Long  Wavelength  Isotropic  Solutions 


For  wavelengths  that  are  long  compared  to  the  size  of  the  scattering 
center  the  calculation  of  acoustic  velocity  and  attenuation  can  be  derived 
ignoring  the  contributions  due  to  scattering.  Even  though  not  all  of  the 
internal  mechanisms  contributing  to  the  dissipation  processes  are  known, 
the  general  theory  can  be  developed  by  grouping  all  attenuation  into  a 
common  source.  This  is  accomplished  mathematically26  by  the  replacement 
of  the  shear  modulus  (M)  and  the  Lame"  constant  (A)  in  the  elastic 
stress-strain  relationship  by  the  first  order  differential  operators 


M  =»  y  +  y 

3t 


A  =  A  +  A' 


In  Eq.  (1),  the  unprimed  terms  denote  the  elastic  and  the  primed  terms 
denote  the  viscous  (or  attenuative)  contributions. 


Consider  an  acoustic  shear  wave*  moving  in  the  positive  (downward)  x- 
direction.  The  displacement  potential  is  given  by 


ir  m  Ae 


i (wt  -  ax) 


(A10) 


where 


u  is  the  circular  frequency 
t  is  the  time 

a  is  the  complex  wavenumber. 

*The  derivation  for  the  dompressional  wave  propagation  is  similar  to  that 
presented  here,  the  difference  being  the  introduction  of  the  (A  +  2M)  op¬ 
erator  in  the  wave  equation,  Eq.  (All),  rather  than  the  shear  modulus  (M) 
operator.  The  displacement  potential  must,  of  course,  be  associated 
with  the  compressional  rather  than  the  transverse  wave. 
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This  displacement  potential  obeys  the  wave  equation 


P 


if i 

3t2 


M  V2  , 


(All) 


where  p  is  the  density  of  the  medium  and  V2  is  the  Laplacian  operator. 
Substitution  of  Eqs.  (A9)  and  (A10)  into  Eq.  (All)  and  carrying  out  the 
operations  gives  the  relation 

pou2  =  a2  (y  +iunO  .  (A12) 


For  an  attenuative  medium,  the  complex  wavenumber  can  be  written  as 

a  =  k  -  it  ,  (A13) 

where 

k  is  the  running  vector,  ai/c 
c  is  the  wave  velocity 
x  is  the  generalized  attenuation. 


Making  this  substitution  into  Eq.  (A12)  and  separating  the  real  and 
imaginary  components  yields  the  simultaneous  equations 


pu2  «=  y(k2  -  t2)  +  2y'u)kt 
0  =  ojy'(k2  -x2)  -  2ykx  . 


(A14) 


Setting 


these  equations  have  the  solutions 
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T 


+  R*u>2 


k2 


•L  p^2  [^  ♦  r2o)2  *  i 

2  y  l  1  +  R2u2 


(A15) 


where  now  both  y  and  R  may  be  frequency  dependent.  The  procedure  from 
this  point  is  to  introduce  a  mathematical  model  of  the  medium  composed  of 
elastic  and  viscous  constants  that  are  frequency  independent,  and  then  fit 
the  results  to  experimental  data. 

Maxwell27  suggested  that  viscoelastic  materials  could  be  represented 
by  an  elastic  element  (spring)  in  series  with  a  viscous  element  (dashpot) 
as  shown  in  Figure  A3.  Voigt28  placed  these  elements  in  parallel  as  in 
Figure  A4.  Detailed  mathematical  analyses  of  these  as  well  as  more  com¬ 
plex  models  are  covered  in  standard  texts29-31  and  will  not  be  discussed 
here.  However,  it  can  be  shown  that  for  the  Maxwell  model 


where 
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,  and 

4- 

i. 


i 


Figure  A3.  Maxwell's  mechanical 
representation  of 
viscoelastic  solids. 


Figure  A4.  Voigt's  mechanical 
representation  of 
viscoelastic  solids. 


For  the  Voigt  model 


T  =* 


I; 

[ 

5 


k2  - 


where 

y  -  Ey 

Ru>  ■  (Ey/nyju  . 
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Very  few  solids  behave  like  either  the  Maxwell  or  the  Voigt  model. 

However,  since  more  complicated  models  become  extremely  involved  mathe¬ 
matically,  and  because  models  specifying  a  greater  number  of  parameters 
require  more  experimental  evidence  substantiate  their  validity,  using 
a  single  Maxwell  or  Voigt  element  is  a  convenient  method  of  obtaining  a 
first  order  approximation  of  a  viscoelastic  solid's  mechanical  properties. 


B.  Short  Wavelength  Solution 

To  estimate  the  contribution  due  to  scattering,  consider  the  problem 
shown  in  Figure  AS.  An  acoustic  plane  wave  propagating  in  an  isotropic 
elastic  solid  encounters  a  spherical  scatterer  which,  for  this  discussion, 
is  also  assumed  to  be  elastic.  Both  media  are  characterized  by  their 


INCIDENT 


SPHERICAL 

SCATTERER 


SCATTERED 

WAVE 


Figure  AS.  Representation  of  the  problem. 


densities  (p^),  shear  moduli  (yj)  and  Lame'  constants  (X^) .  Elasticity 
theory  defines  the  velocities  as 


(A18) 


where  u>  is  the  circular  frequency  and  and  are  the  wavenumbers  of 
the  congressional  and  transverse  wave,  respectively.  The  displacement 
can  be  written  in  terms  of  the  displacement  potentials  $  and  $  as 

u  ■  grad  <J>  +  cur£  $  .  (A19) 


For  spherical  coordinates, 

u  ■  -7<P  ♦  Vx7x(r\p)  . 


(A20) 


The  displacement  potentials  are  of  the  form 


E  PmCcose)  eiWt  *  (A2D 

mo 


where 


fB(pr)  ■  the  spherical  Bessel  function  of  order  m 

p  ■  either  k  or  K,  depending  upon  the  wave  considered 


PB(cos8)  *  the  Legendre  polynomial  of  order  m. 


When  the  incident  plane  wave  strikes  the  scattering  center ,  some  of  the 
energy  is  scattered  away  and  some  is  transmitted  through  the  scatterer. 
Since  all  the  waves  must  obey  the  same  wave  equations,  the  wave  potentials, 
neglecting  the  time  dependence,  may  be  written  as 


*  Incident  =  £  I  3m  (V )  Pm  (cos6) 

m»o  (A22) 

00 

ip  Incident  =  ^  I  j  (K,r)  P  (cos0) 
sm  m  *  m 

m*=o 


4>  Scattered  =  £  A  h  (k.r)  P  (cosG) 
m  m  1  m 

m=o 

(A23) 

00 

ip  Scattered  *  (K^r)  (cos6)  » 

m«o 


<j>  Transmitted  *  £  Cm  ^m  ^k2r^  Pm  (cos0) 
m*o 

(A24) 

00 

Transmitted  =  ^  D  j  (K~r)  P  (cos0) 

m  m  i  m 

m*o 

where  r  and  0  are  spherical  coordinates.  The  spherical  Bessel  functions  of 
the  first  kind,  jm(kr),  and  third  kind,  h^kr),  are  used  to  assure  that  the 
scattered  wave  is  propagating  away  from  the  scatterer  located  at  the  origin 
of  the  coordinate  system.  In  order  for  the  incident  wave  to  be  plane  and  of 
unit  amplitude,  Icm  and  Ism  must  be  of  the  form* * 


I 


m 


(2m  +  1) 


I 


(A2S) 


where  p  is  k  and  K,  respectively.  The  coefficients  A„,,  Bg,  C,,,  and  Dm  are 
determined  by  the  boundary  conditions  at  the  scattering  surface.  For  a 
solid-solid  interface  assuming  strong  coupling,  these  boundary  conditions 
are 


(a)  Radial  displacement  continuous 

(b)  Tangential  displacement  continuous 

(c)  Radial  stress  continuous 

(d)  Tangential  stress  vanishes. 


At  the  boundary,  r  ■  a,  these  conditions  can  be  written  mathematically  as 


(a) 

ur(inc) 

+  ur(scatt) 

ur(trans) 

(b) 

u6(inc) 

+  u6(scatt) 

u0(trans) 

(O 

°rr(inc) 

+  0rr(scatt) 

3  (7  - 

rr(trans) 

(d) 

°r6(inc) 

+  °r0(scatt) 

ar0(trans)  . 

(A26) 


In  terms  of  the  displacement  potentials,  the  stress  and  displacement  com¬ 
ponents  of  Eq.  (A25)  are 


o _  «  pu> 

rr 


r 

where 


jr «•!?(?»  *)]| 


Q  "  5lb!?(sin0  w)  * 


(A2S) 


jr 
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The  scattering  cross  section  is  determined  by  calculating  the  ratio  of 
the  time  rate  at  which  energy  is  scattered  out  of  the  wave  by  an  obstacle  of 
radius  "a"  to  the  total  energy  of  the  incident  wave.  The  scattered  energy  is 
equal  to  the  energy  being  carried  away  by  the  scattered  wave  across  a  spheri¬ 
cal  surface  of  radius  "b"  >  "a".  Love*3  gives  this  scattered  energy  as 


^scattered  ■// 
A 


xr 


3u 

x 

w 


♦  a 


yr 


3t  +a 


3u 


_ z 

zr  3t 


scattered 


dA  , 


(A29) 


where  axr,  a^,  and  ozr  are  the  stress  components  acting  in  the  three  rectan¬ 
gular  axes  on  a  surface  normal  to  the  radius  vector  r,  and  the  u's  are  dis¬ 
placements.  Both  the  o's  and  the  u's  are  usually  complex,  so  care  must  be 
exercised  to  assure  that  the  final  expression  for  the  scattered  energy  is 
real.  Assuming  a  time  dependence  of  eiu,t  for  the  displacement  potentials 
and  using  the  spherical  symmetry  of  the  scattered  wave,  Eq.  (A29)  can  be 
rewritten  in  spherical  coordinates  as 


scattered 


ff  { |°rr  ur  *  °9r  V°zrul| 


r*b 


(ASO) 


-  [°rrar  *VU0  M’xrui  »i«“8 

l  ’  ’ scattered 


scattered 
wave 


where  all  the  terms  are  defined  by  Eq.  (A26)  and  the  asterisk  denotes  the 
complex  conjugate.  The  integral  in  this  equation  can  be  evaluated  by  the 
substitution  of  Eqs.  (A23)  and  (A24)  and  using  the  identity 


/*CBf)  g  sin0d6  - /*f (fig)  sinSde  -  -/*  ||  |f  sin0d6 


(A51) 
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to  determine  the  recursion  relations  and  the  orthogonality  of  the  Bessel 
functions  and  the  Legendre  polynomials.  Equation  (A 31)  is  valid  for  any 
two  arbitrary  functions  f  and  g  and  the  energy  lost  from  the  incident  wave 
due  to  scattering  from  the  spherical  scatterer  is 


(A32) 


where  Ag,  and  Bg,  are  the  magnitudes  of  the  scattered  wave  potentials  and 
must  be  evaluated  at  the  scattering  surface  by  the  boundary  conditions. 


The  total  energy  of  the  wave  is  determined  by  calculating  the  energy 
in  the  incident  plane  wave  that  is  being  transported  through  a  unit  area 
normal  to  the  propagation  direction.  The  displacement  of  a  plane  wave 
of  unit  amplitude  traveling  in  the  positive  x-direction  can  be  written  as 

t  -  z  ei(wt  *  Kl2)  .  (A33) 

From  elasticity  theory,  the  stress  can  be  calculated  from  the  relation 


(A34) 


1  i-j 

0  if»j 


For  the  wave  described,  ojj  ■  0  for  all  ij  except 


o 


zz 


R1 


(A35) 


P,w* 


The  energy  flux  through  a  unit  area  is  then  given  by 

E  ■  i  |(o  u*  ♦  a  u*  ♦  a  u#)  -  (o  *u  ♦  a  *u  ♦  o  *u  )  \  ■  ~— 
inc  2  \v  xz  x  yz  y  zz  z'  v  xz  x  yz  y  zz  z'J  Kj 

The  scattering  cross  section  is  defined  as 


(A36) 


"scattered 

Etotal 


(A37) 


Therefore,  for  a  transverse  plane  wave  encountering  a  liquid-filled  spherical 
scatterer 


''shear  ■  4"  E 


STT  ET  l». 

m*o  i  1 


•|!]  • 


(A38) 


A  similar  calculation  for  the  problem  of  a  plane  congressional  wave  striking 
the  liquid-filled  cavity  yields  the  similar  expression 


^comp"  **  ^  2m  ♦  1 

r  m«o 


Am|2  ♦  mC**l) 


Jr 


(A39) 


The  scattering  cross  section  times  the  square  of  the  frequency  (fv2)  as  a 
function  of  the  circumference  to  wavelength  of  the  scatterer  has  been 
solved  using  the  computer  program  SCTTER  given  in  Appendix  B  and  the  results 
shown  in  Figures  13-15  of  the  text. 
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APPENDIX  B 


COMPUTER  PROGRAMS 


PROGRAM  SOLID 


ATTACH (GRAFIX) 

LIBRARY (GRAF IX) 

FORTRAN. 

LGO(LC>3000) 

PROGRAM  SOLIOWFC INPUT, OUTPUT. TAPE5=INPUT,TAPE6=OUTPUT,T APE99) 
IMPLICIT  COMPLEX  (A«F*G)*REAL(K) 

COMPLEX  BET A 

DIMENSION  0(4*5) •  A(4,5),  F(4),  KK(4),  ENG(8),  CP(8,1000) 

COMMON  N,  NN,  G.  F 

CALL  PLOTS 

CALL  PLOT  ( 0  •  t2«  .  -*3) 

NNX*Q 

READ (5* 100) NNN. NP 

100  FORMAT  (212) 

101  READ(5,102)C0,CR,CL,CT,RHC0,RHC,CHI0,DCHI,CHIF 

102  F0RMAT(4F5.Q,2F7.4,3F8.4) 

XI=DCHI/5. 

XJ-OCH I*  8 • 

PI=3. 141592654 

RsRHO/RHOO 

K0=1./C0 

KR*1./CR 

KL*1./CL 

KT*1./CT 

LL*IFIX( (CHIF-CHIO) /OCHI*. 01) 

CHI0*CHI0*PI/180 . 

OCH I=OCH I*PI/ 18  0 • 

CHIF=CHIF*PI/180. 

CHI*CHI0 
00117IL-1 *LL 
CHEK1=0. 

CMEK2=  0 • 

IF (CHI) 104*105*104 

104  IF (I L- 1)106*106*105 

105  CHI=CHI*OCHI 

106  BETA=KO*SIN(CHI) 

AO*KO*COS (CHI ) 

AR*CONJG(CSQRT(KR**2-8ETA**2)  ) 

AL*CONJG (CSQhT (KL**2-BETA**2) ) 

AT *CONJG (CSQRT (K  T**  2-BET  A** 2) ) 

N»4 

NN*N+1 

00  107  I« 1 , N 
DO  107  M=t»  NN 

107  G(I*M)sQ. 

G(l,l)sl. 

6(1,2)=- BET A/ AO 
G(l, 3) sAL/AO 
G(l»4)*-G(l,2) 

G(l*5) =1. 

G(2*l)=l. 

G(2,2) sAR/BETA 
G(2,3) *-l. 

G  (2, 4)* AT/BET A 
G (2,5) *-l. 

G(3tl)*-(1.-2.M8ETA**2/KR**2> ) 

G«3,2)*2.*8ETA*AR/KR**2 


G(3,3)*(i.-2.*(BETA**2/KM*2>  >*R 

G(3,4)*2.*R*8ETA*AT/KT**2 

G(3,5)*-G(3,l) 

G(4,1)*2.*A0*BETA/KR**2 

G(4,2>=-G(3,l> 

G(4»3)*AL*G(3*4)/AT 

G(4,4»=-G(3,3> 

G(4,5)=G(4,1> 

CALL  SOLN 
DO  108  I J=1 « 4 
108  KK(IJ)=CAaS(F(IJ> 1**2 
ENG(1)=KK(1) 

ENG(2)=REAL( (AR/A0I*KK(2) ) 
ENG(3)=REAL(RMAL/A01  *KK(  3)  ) 
ENG(4)=REAL(R*< AT/AO) *KK (4) ) 

00  200  I  J  =  1  ,  4 
CHEK1=CHEK1*ENG ( I J) 


200  CONTINUE 

00  109  1*1, N 
00  109  H*1,NN 


109  G (I 

,M> 

*0. 

G(1 

*1) 

=AO/BETA 

GC1 

,21 

=  -l. 

G  ( 1 

,3) 

=  AL/BET  A 

G  (1 

,41 

*1. 

Gil 

,5> 

=  1. 

G  (2 

,1) 

*BETA/AR 

G  (2 

,2) 

=  1. 

G  (2 

,3) 

=-G (2, 1) 

G  (2 

,4) 

= AT/AR 

6(2 

,51 

*1. 

G  (3 

,1) 

*-(l.-2.*(BETA**2/KR* 

G  (3 

,2) 

*2.*BETA*AR/KR**2 

G  ( 3 

.3) 

=  RM1.-2.MBETA**2/KT 

G  (3 

,4) 

*2.*R*8ETA*AT/KT4#2 

G  (3 

,5) 

=G(3,2) 

G  (4 

.1) 

*2.*A0*8ETA/KR**2 

G  (4 

,2) 

*-G(3,l> 

G  (4 

,3) 

*AL*G ( 3,4) /AT 

G  (4 

,4) 

=-G (3 , 3) 

G  (4 

,5» 

*G (3, 11 

CALL  SOLN 
00  110  IJ*1,4 

110  KK(IJ)*CABS(F  C I J ) 1**2 
ENG(5)=REAL((A0/Afi)*KK(l)) 

ENG (6) *KK(2) 

ENG (7) *REAL(R*(AL/AR)*KK(3) ) 

ENG(8)*REAL(R*(AT/ARI  *KK(41I 
00  201  I J=5 , 8 
CHEK2*CHEK2  ♦  ENGCIJ) 

201  CONTINUE 

THETA*CHI  *  180. /PI 
IF(IL-1)111,111,U5 

111  MRITE(6, 112) 

112  FORMAT ( 1H1I 

WRITE  (6,1131  CO,  CR,  CL,  CT,  RHO,  RHOO 

113  FORMAT  ( 3X«  5HC0  *  ,  F5.0,  3X ,  5HCR  *  ,  F5.0,  3X,  5MCL  *  , 

t  FS.O,  3X,  5MCT  *  ,  F5.0,  3X,  6HRH0  *  ,  F4.2,  3X,  7HRM0I  a  , 
S  F4.2,  /» 
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o  o 


WRITE  (6,  114) 

114  FORMAT  (8X,  33H - COMPRESSIONAL  WAVE  INCIDENT---,  6X, 

S33H - TRANSVERSE  HAVE  INCICENT - ,  /,  IX,  5HTHETA, 

S3X,  2HA2,  5X ,  2HB2,  5X,  2HA3,  5X,  2H83,  4X,  5HCHEK1,  EX, 

I2HA2,  5X ,  2HB2 ,  5X,  2HA3,  5X,  2H83,  4X,  5HCHEK2) 

A3  AND  B3  ARE  THE  TRANSMITTED  COMPRESSIONAL  ANO  TRANSVERSE 
HAVE  ENERGIES  RESPECTIVELY. 

115  WRITE  (6,116)  THETA,  (ENG(Nf),  NM=1,4),  CHEK1,  (ENG(NM),  NM=5,5), 
SCHEK2 

116  FORMAT  (IX,  F5.2,  4(2X,  F5.3),  2X,  F5.3,  4X,  4(2X,  F5.3) ,  2X, 
JF5.3I 

00  10  IP  =  1 » 6 

CP( IP, IL )  =  ENG ( IP )  *  10. 

10  CONTINUE 

117  CONTINUE 
IP  *  0 

DO  55  J= 1 , 2 
YP=6. 

00  45  I Q  =  1 , 4 

RI*IQ 

IR=IQ-1 

IF (IR.EQ.3)  IR=5 
IP  *  IP»1 
XN  x  XI 

CALL  PL0T(XN,CP(IP,1) ,3) 

00  35  IL  =  2 , LL 
XN  a  X  N  +  X I 

CALL  PL0T(XN,CP(IF,IL)  ,2) 

35  CONTINUE 

XN  a  X J/4.*RI 
ISaIQ»10 
YP=  YP- .25 
TS*YP*.07 

CALL  SYMBOL (2.,YS,.14,IR,0.*-1) 

CALL  SYMBOL (2.17,YP,.14,1H-,0.,1) 

00  40  ILSIS,LL ,40 

CALL  SYMBOL(XN,CP(IP,IL) ,  .14, IR,  0. ,-l) 

XN  a  XN+X J 
40  CONTINUE 
45  CONTINUE 

CALL  PLOT  (18.,  0.,  3) 

CALL  PLOT  (0.,  0.,  2) 

CALL  PLOT  (0.,  10.,  2) 

CALL  PLOT(18.,10.,2) 

CALL  PLOT  (18.,  0.,  2) 

XP*16. 

YP»  0.1 
Y2*0. 

00  451  Lai, 8 

CALL  PLOT  (XP,  YP,  3) 

CALL  PLOT  (XP,  Y2,  2) 

XPaXP  -  2. 

451  CONTINUE 
YP«1. 

XPaQ.l 

X2«0. 

OC  452  Lai, 9 

CALL  PLOT  (XP,  YP,  3) 
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CALL  PLOT  (XZ,  YP,  21 
YP=YP*1. 

452  CONTINUE 
YP=-0. 07 
XP=-.42 
YN*0. 

00  453  L=l, 11 

CALL  NUMBER  (XP, YF,  .  14, YN ,  0 . , 1 ) 

YN=YN*.l 

YPsYP+1. 

453  CONTINUE 

CALL  SYMBOL  (-. 65 ,4 . 1 , .14 , 15HREL AT  IVE  ENERGY, 90. ,15) 

IP  (J  .EG.  1)  X  P=  6 . 3  6 
IF  (J  .EQ.  2)  XP=6 .54 
IF ( NP. EQ • 1 )  GO  TO  230 

CALL  SYMBOL  ( 1 . 7 , 10 .4 , .21 , 8 1HREL AT IVE  ENERGY  DISTRIBUTION  OF  ACOUS 
STIC  NAVE  INC IOENT  ON  PERM AFRO S T- ICE  BOUNDARY, 0. , 61) 

IF  (NNX  .EQ.  0)  CALL  SYMBOL  (XP+.06,  10.10,  .14, 

S 15H0TT  AHA  SAND-I CE  ,  0  .  ,  15 ) 

IF  (NNX  .EQ.  1)  CALL  SYMBOL  ( XP , 1 0. 1 0 ,. 14 , 16HHAN0VER  SILT-ICE, 

S  0 . • 16) 

IF  (NNX  .EQ.  2)  CALL  SYMBCL( XP-  .06, 10.1, . 14, 17HG00DRICH  CLAY-ICE, 
SO. ,17) 

GO  TO  240 
230  CONTINUE 

CALL  SYMBOL(0.8,10.4,.21,91HRELATIVE  ENERGY  DISTRIBUTION  OF  ACCUST 
SIC  HAVE  INCIDENT  CN  BOUNDARY  OF  TWO  PERMAFROST  L A YERS , 0 . , 9 1 ) 

CALL  SYMBOL (XP-. 54, 10.1, • 14*25H0TTAWA  SAND-GOCDRICH  CLAY, 0. ,25) 

240  CONTINUE 

IF  (J  .EC.  2)  GC  TO  454 

CALL  SYMBOL (999. ,999. , .14, 3 OH  (CCMPRESSICNAL  HAVE  INCIDENT) , 0.  ,30) 
GO  TO  455 

454  CONTINUE 

CALL  SYMBOL  ( 999. , 999 . , . 14 , 27H  (TRANSVERSE  HAVE  INCIDENT) , 0.  ,27) 

455  CONTINUE 

CALL  SYMBOL (2.35,5.75, .14, 2 3HREFLECTEC  CCMPRESSICNAL, 0 . ,23) 

CALL  SYMBOL  ( 2 . 35, 5 .5 , .14 , 15HREFLECTED  SHEAR, 0., 15) 

CALL  SYMBOL  ( 2. 35 ,5 .25 , . 14 ,25H TRANSM I TTED  COMPRESSION AL , 0 ., 25) 

CALL  SYMBOL  (2 . 35 ,5 . 0 , .14 , 1 7HT RANSMIT TED  SHEAR, 0. ,17) 

XN*10. 

XP*1.9 

YP*-0.24 

CALL  NUMBER  (-. 04, YP, .14, 0. ,0 . ,-l) 

00  456  L=l,9 

CALL  NUMBER  ( XP, YF, . 14 , XN , 0 . , - 1 ) 

XN*XN* 10 . 

XP»XP42. 

456  CONTINUE 

CALL  SYMBOL  ( 7. 5 8,-0 .48 , . 14 ,2 4HINCIDENT  ANGLE  (DEGREES),  0.,24) 
CALL  PLOT (24. ,0 , »-3) 

55  CONTINUE 

NNX  ■  NNX  ♦  1 

IF  (NNX  -  NNN)  101,  118,  118 
118  CONTINUE 
X*0. 

Y«t. 

00  130  L*l»4 
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LP=L 

IF (LP.  EQ  .4)  LP  =  5 
00  120  LX-1,11 

CALL  SYM80L<X,Y, . 14 ,LP , 0  .  , - 11 
Y=Y*1. 

120  CONTINUE 
X=X*1. 

y=o. 

130  CONTINUE 

CALL  PLOTCO.,0.,9991 

STOP 

END 

SUBROUTINE  SOLN 

IMPLICIT  COMPLEX  (A,  F,  G),  REAL  <K> 

COMPLEX  SUMPR 

DIMENSION  G  (  4 , 5 1  ,  A(4,5>,  F  C  4 ) 

COMMON  N,  NN •  G,  F 

C  THE  auxillary  matrix  for  the  solution  of  the  matrix  by  the  crout 

C  METHOD  FOLLOWS. 

DO  104  1  =  1,  N 

104  A(I* II =G (I, II 

00  111  1=2, N 

00  111  M=2,NN 

A(1,M)=G(1,M)/A(1,1) 

SUMPR=CMFLX ( 0 • , 0  .  I 

IF ( I-Nl  105,  107,  109 

105  11=1-1 

00  106  MM= 1,11 

106  SUM PR=  SUMPR*  A (I, MM) *  A (MM ,  M 1 
A 1 1, Ml =  <G d,MI-SUMPRI/A<I,II 
GO  TO  111 

107  11=1-1 

00  108  IJ=1,II 

108  SUMPR=SUMPR*A (I,IJ)*A(IJ,I) 

A(I*II=G(I, II- SUMPR 

GO  TO  111 

109  MM=M-1 

00  110  I J=1 » MM 

110  SUMPR=SUMPR* A (I,IJ)*A(IJ,M) 

ACI,M)=G(I, Ml -SUMPR 

111  CONTINUE 

C  THE  FOLLOWING  DETERMINES  THE  FINAL  MATRIX  BY  THE  CROUT  METHOO- 

C  FOR  THE  FROBLEM  CONSIDERED. 

F ( N! =A ( N , NNI 
MM=N-1 

00  115  M=1 , MM 

L*N-M 

ML=L*l 

SUMPR=CMFLX ( 0. ,0.1 

112  IF (N-MLl 114,  113,  113 

113  SUMPR=SUMPR*A<L,MLI*F<MLI 
ML=ML»1 

GO  TO  112 

114  F (LI =A (L , NNI -SUMPR 

115  CONTINUE 
RETURN 
ENO 
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PROGRAM  TRNTION 


ATTACH (GRAFIX) 

LIBRARY (GRAF  IX) 

FORTRAN. 

LGO(LC*6000) 

PROGRAM  TRNT ICN( INPUT . OUTPUT, TAPE5= INPUT .TAPE6=0UTPUT) 
THIS  PROGRAM  CALCULATES  THE  REFLECTION  COEFFICIENT  FROM  A 
TRANSITION  LAYER  BETWEEN  THC  ELASTIC  MEOIA. 

IMPLICIT  REAL (K,  N) 

COMPLEX  I ,GM1 *GP1 * A2 

DIMENSION  HOLD (2) ,X (500) , YC (500) ,YT(500) 

CALL  PLOTS 
CALL  PLOT (0..1..-3) 

116  CONTINUE 

100  READ (5*101)CQ,CR,CC,CS,RHOO,RH02,N,MNX 

101  FORMAT (4F5.0.2F5. 3, F3. 0,13) 

IF (EOF ,  5 )  300.102 

102  CONTINUE 
IF(HNX.EQ.O)  GO  TC  1 
PS=40. 

PL*  .025 
GO  TO  2 

1  CONTINUE 
PS=20. 

PL*. 05 

2  CONTINUE 
NT  =  0  . 

I=(0.,1.) 

K0H*0« 

00  114  L=1.500 
K0H=K0H+ • 02 
NT*NT+1. 

104  KRH=K0H*C0/CR 
KCH*KOH*CO/CC 
KSH*K0H*C0/CS 
IF(N)105,105,106 

105  ALPH=0 . 

ALPT*0. 

GO  TO  107 

106  XN=1 ./N 
ALPH=(CC/C0) **XN-1. 

ALPT*( CS/CR) **XN-1. 

C  KLH  ANO  KTH  OEFINE  THE  TRANSITION 

107  KLH=K0H/ (l.+ALPH) **N 
KTH*KRH/(l.tALPT)**N 
NM1=N- 1 . 

C  FIRST  CALCULATE  FCR  COMPRESSI ONAL  INCIDENT 
P*(1.-(N-1.)*ALPH)/(1.*ALPH)**NM1 
Q*(1.-(N-1.) *ALPH)**2 

R*(2.*ALPH*N/K0H) M1.-0.5*  (N-1.)*ALPH)*( (1 . ♦ALPH ) **N ) 

XI*  (RHO2/RHO0) * (KOH/KCHI*P-Q+R*TAN (KLH) 
X2*RMQ-(K0H/KCH)»P) *T  AN ( KLH) 

X3*Q-R*TAN (KLH) 

X4* (KOH/KCH) *P*TAN ( KLH) 

Y1*(RH02/RHOO)*(KOH/KCH) *P*C-R*TAN (KLH) 

Y2*R*(Q* (KOH/KCH )*P)*T AN (KLH) 

GM1* (X1*I*X2)/(X3*I*X4) 

GP1*(Y1*I*Y2)/(X3«I*X4) 
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A2*(GN1-2.*I*N*ALPH/KOH>/(GP1-2.*I*N*ALPH/KOH> 

HOLD (1)=CABS(A2) 

Z*(RH02*CC)/ (RHOO*CO) 

TEST*(l.-Z)/(i.»Z) 

IF(HOLD(1)-TEST)201,201,200 
2U0  HOLD ( 1 ) =  TEST 

201  CONTINUE 
ENGC=H0LC(1)**2 

C  NOW  CALCULATE  FOR  TRANSVERSE  HAVE  INCIDENT 
P=I1.-(N-1.) *ALPTI/(1 .♦ALPT)**NM1 
Q«U.-(N-1.I*ALPT»**2 

R«C2.*ALPT*N/KRH)*(l.-0.5MN-l.)*ALPTIM(i.*ALPT>»*N) 
Xl=(RHO2/RHO0>* (KRH/KSH) *P-Q*R*TAN (KTH) 

X2*R* (Q-( KRH/KSH ) *P ) *T AN (KTH) 

X3=Q-R*TAN(KTH) 

X4* (KRH/KSH) *P*T  AN(KTH) 

Yl* (RH02/RHOO) *( KRH/KSH) *P*Q-R*TAN(KTH) 

V2=R*(QMKRH/KSH)*P)*TAN(KTHI 

GHI*(X1»I*X2»/(X3*I*X4) 

GP1*CV1»I*Y2>/(X3*I*X4) 

A2= (GM1-2. * I*N*ALFT/KRH)/(GP1-2.*I#N*  ALPT/KRHI 
HOLD (2)sCABS(A2) 

Z=(RH02*CS)/ (RHOO»CRJ 
TEST=(1.-Z)/(1.*Z> 

IF (HOLD (2 1 -TEST)  203,203,202 

202  HOLD(2) =  TEST 

203  CONTINUE 
ENGT=HOLOI2)**2 
IFINT-l.) 108,108,115 

108  CONTINUE 

IF(HNX.EQ.O)  WRITE (6,111) 

IF(HNX.EC.l)  WRITE (6, 121) 

IF ( HNX • EQ • 2 )  WRITE (6,131) 

111  FORHAT (1H1,16X,15H0TTAWA  SAND-ICE) 

121  FORMAT (1H1,15X,16HHANCVER  SILT-ICE) 

131  FORMAT (1H1,15X,17HG000RICH  CLAY-ICE) 

HRITE(6, 109)  N 

109  FORMAT (10X,27HTRANSITICN  OF  FORK  (1*AX>** 

1 ,F2.0,/,12X,17H. .COMPRESS ICNAL. . ,3X, 17H. .. TRANS VERSE  .... ) 

MRITE(6, 112) 

112  FORMAT  (5X,3HK0H*4X,7HCMP  AMP,3X,7HCMP  ENG,3X,7HSHR  AMP,3X,7HSHK  EN 
1G) 

115  NRITE(6,113)KOH,HCLO(1),ENGC,HOLD(2) ,ENGT 

113  FORHAT  (1X,F8.3,4(3X,F7.4)  ) 

X(L)*K0H 

YC(L)=ENGC*PS 

YT(L)=ENGT*PS 

114  CONTINUE 

CALL  PLOT (X ( 1) ,YC(1> ,3) 

00  10  J*2,5QQ 

CALL  PLOT(X(J),YC(J),2) 

10  CONTINUE 

DO  15  J=20,500,20 

CALL  SYMBOL (X(J) , YC ( J ) , . 1 36 ,0 , 0 . , - 1) 

15  CONTINUE 

CALL  PLOT(X(l),YT(l) ,3) 

00  27  J»2,5Q0 


CALL  PLOT (X(J) » Y  T (Jl ,21 
CONTINUE 

00  46  J3 10, 49 0,20 

CALL  SYMBOL (X  ( J) , YT ( J) ,. 136 ,1 , 0 .  ,-l) 

CONTINUE 

CALL  PLOT (10. ,0 . ,3) 

CALL  PLOT ( 0.  ,0. *2) 

CALL  PLOT  (0. ,8. *2) 

CALL  PLOT (10*  ,8. ,2) 

CALL  PLOT (10*  ,0. ,2) 

XP*9. 

YP*0.1 

VZ*0. 

00  48  J  =  1 ,9 
CALL  PLOT  (XP*YP«  3) 

CALL  PLOT (XP,YZ«2I 
XP=XP-1 . 

CONTINUE 
YP=  1 . 

XP=0.1 

XZ*0. 

00  60  J  =  1 , 7 
CALL  FLCT(XP,YP,3) 

CALL  PLOT  ( XZ  *  YP, 2) 

YP=VP»1. 

60  CONTINUE 
YP=-0.07 
XP*-.7 
YN=0. 

00  65  J  =  1 «  9 

CALL  NUMeER(XP,YP,.14,YN,0.,3) 

YN=YN*PL 

YPsYP^l. 

CONTINUE 

CALL  SYM80L(-.89,3.1, .14, 15HRELATIVE  ENERGY, 90 ., 15) 

CALL  SYMBOL (2,66, 8,40, *21,2 5HTRA NS ITION  CF  FORM  (1*AX> ,0. ,25) 

CALL  NUMBER (999. ,8. 56 , .14 ,N , 0 .  ,-il 

IF (HNX.EQ.O)  CALL  SYMBOL  I  4. 10 , 8 . 1 ,  .14 , 15H0TTAMA  SANO-ICE  ,0.,15) 
IF (HNX.EQ.l)  CALL  SYMBOL ( 4. 04 , 8 . 1 , .14 ,16HHAN0VER  SILT-ICE  ,0.,16) 
IF (MNX.EQ.2)  CALL  SYMBOL (3. 98, 8.1, .14, 17HG000RICH  CLAY-ICE ,0 ., 17) 
CALL  SYMBOL ( 6.0, 7.0,. 136, 0,0. ,-l) 

CALL  SYMBOL <6. 17, 6. 93,. 14, lH-,0., 1) 

CALL  SYMBOL (6.35,6.93,. 14 ,22H COMPRESS IONAL  INCIDENT, 0. ,22) 

CALL  SYMBOL (6. 0,6. 75,. 136, 1, 0.  ,-l) 

CALL  SYMBOL (6. 17,6.68,. 14, 1H-,0.,1) 

CALL  SYMB0L(6.35,6.68,.14,14H SHEAR  INCIDENT, 0. ,14) 

XN*0. 

XP*-.04 
VP*-0.24 
00  75  J*l,10 

CALL  NUMBER (XP, YP , .14,XN,0.,-1) 

XN«XN+1. 

XP«XP»1. 

75  CONTINUE 

CALL  NUMBER (9.9, YP, .14,XN,0.,-1> 

CALL  SYMBOL (4.82«-.48,.14*lHK,0.,l) 

CALL  SYMBOL (4. 944, -.5 36,. 112, 1H0,0.,1) 
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CALL  SYMBOL (5. 064, -.48,.  14,1HH,0.,1> 
CALL  PLOT (14. , 0.  »-3) 

60  TO  116 
300  CONTINUE 

CALL  PLOT ( 0 . «  0. , 9991 

STOP 

ENO 


PROGRAM  SCATER 


FORTRAN. 

L60* 

PROGRAM  SCATER C INPUT, OUTPUT, T APE5=INPUT, TAPE6*0UTPUT) 

IMPLICIT  REAL  (J,KI 

COMMON  KL1,KS1,KL2,KS2,RH01 ,RH02, SIM, CIM, GAMMA. AK,J (400) • NMAX . 

1 Y (400) ,CL1,CS1,CL2,CS2 

THIS  PROGRAM  CALCULATES  THE  SCATTERING  CROSS  SECTION  TIMES  THE 
FREQUENCY  SQUARED  FOR  A  SHEAR  AND  COMPRESS  ION AL  PLANE  HAVE 
IMPINGING  ON  AN  ELASTIC  SPHERE  EHBEOOEO  IN  AN  ELASTIC  MEDIUM. 

THE  VALUE  IS  CALCULATED  AS  A  FUNCTION  OF  THE  OIMENSIONLESS 
PARAMETER  KA  WHERE  K=2.*PI*FREQUENCY/VELCCITY  AND  A  IS  THE  SPHERE 
RAOIUS .  THE  CALCULATION  IS  MADE  IN  THE  C.G.S.  SYSTEM.  THE 
SCATTERING  CROSS  SECTION  HAS  UNITS  OF  DB-CM**2.  CS1  AND  CS2  ARE 
THE  SHEAR  ANQ  CL1  AND  CL 2  ARE  THE  COMPRESSIONAL  VELOCITIES  IN  THE 
MEOIUM  AND  SCATTERER  RESPECTIVELY.  RHOl  AND  RH02  ARE  THE 
RESPECTIVE  DENSITIES. 

MRITE<6.97) 

97  FORMAT ( 1H1 , 45X  « 12HHANCVER  SILT) 
REA0(5,98)CL1,CS1,CL2,CS2,RH01,RH02 

98  FORMAT (4F8.0.2F6.3) 

CL1*CL1*10C. 

CS1*CS1M00. 

CL2*CL2*100. 

CS2»CS2*100. 

MRITE(6,99) 

99  F0RMAT(5X,2HKA,7X,11HGAMMA, SHEAR, 7X, 10HGAMMA.COMP) 

100  READ (5 • 101) KA 

101  FORMAT (F7. 2) 

IF (EOF. 9)1 04 .102 

102  CONTINUE 

FIRST  THE  SHEAR  SCATTERING  IS  CALCULATED 
KS1«KA 

KL1*KS1*CS1/CL1 
KS2*KS1*CS1/CS2 
KL2*KS1*CS1/CL2 
SIM*CSl 

IF  CIM«0.,  THEN  A  SHEAR  WAVE  IS  INCIOENT. 

CIM«Q. 

CALL  SOLN 
GAMMAS«GAMMA 

GAMMAS  IS  THE  SHEAR  SCATTERING  CROSS  SECTION  *  FREQ»'2. 

NON  CALCULATE  THE  COMPRESSIONAL  SCATTERING  CROSS  SECTION. 

KL1*KA 

KS1«KL1*CL1/CS1 
KL2>KL1*CL1/CL2 
KS2*KL1*CL1/CS2 
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C  IF  SIMaO.,  THEN  A  CONPRESSICNAL  HAVE  IS  INCIOENT 
SIH«0. 

CIM«CL1 
CALL  SOLA 
GAHHAC«GAMHA 

C  6AHHAC  IS  THE  COFFRESSlONAL  SCATTERING  CROSS  SECTI0N*FREQ**2 
HRITE(6,103) KA,G4MMAS  «GAMMAC 
103  F0RMAT(2X,F7. 2, 2I5X, £13.4)1 
GO  TO  100 
10*  CONTINUE 
ENO 

SUBROUTINE  SOLN 
IMPLICIT  REAL ( J ,  K) 

COMPLEX  A,F,G,I,$LMPR 

DIMENSION  Jl(400)  ,J2(40U>  , J 3 ( 4 00 > , J4 <400 > , Y1 t 4 00 >  ,Y2(400)  ,G<4,5>  , 
1A  (4 » 5) »F  (4) »  K K ( 4 ) 

COMMON  KLl,KSi,t<L2,KS2,RHCl,RH02, SIM, CIM, GAMMA, AK,J(400> , NMAX , 

1 Y ( 400) ,CL1,CS1,CL2,CS2 

THIS  ROUTINE  SOLVES  THE  BOUNDARY  CONDITIONS  FOR  THE  AMPLITUDES  OF 
THE  SCATTERED  HAVE  HHICH  ARE  USED  TO  CALCULATE  THE  SCATTERING 
CROSS  SECTIONS. 

Is  ( 0  *  ,  1 .  ) 

PI»3. 141 592654 
CALL  UPRLIM 
NM*NMAX+1 
AKaKLl 
CALL  SPHBES 
DO  100  L«1,NM 

100  JKD-JUI 
AKsKSl 
CALL  SPHBES 
00  101  L*1,NM 

101  J2(L)«J(L) 

AK*KL2 
CALL  SPHBES 
DO  102  L»1 , NM 

102  J3(L)sJ(L) 

AK»KS2 
CALL  SPHBES 
00  103  L>1,NM 

103 

NH*NMAX>2 
AK«KL1 
CALL  SPHNEU 
00  111  M«1,NH 

111  Yl(H)aY(M) 

AKaKSl 
CALL  SPHNEU 
DO  112  M*1,NM 

112  Y2(M)»Y (H) 

GAMMA  *0. 

00  127  Maif NMAX 


L 


Nl»4 

N2*N1*1 

00  113  NX=1,N1 
00  113  NY»1 »  N2 
113  G (NX,NY) *0. 

G(1»1)=-(M*J1(M)  -KL1*J1(M*1)) ♦I*(M*Y1(H)-KL1*Y1(N»1) ) 
G(1,2)*(M*(N*1)*J2(M) )-I*  (M*(H*1)*Y2(H)) 

G (1*3>* (N*J3(M) -KL2*J3(M*1) I 
G(i,4)*-N*(M*1)*J4(M) 

G(l,5)sCIM*(M*Jl(M)-KLl*Ji(M+l))-SIN*(M*  (M*l) *J2(M) ) 
G(2,1)*-J1(M>*I*Y1(M) 

G(2,2)=( (M*1)*J2(M)-KS1*J2(M*1)> -I*( (M*1)*Y2(M)-KS1*Y2(N*1)> 

G(2«3) =J3(M) 

G(2,4)=-((M*1)*J4(N)-KS2*J4(M*1) ) 

G(2,5)=CIM*J1(M) -SIM* I (M*1)*J2(M)-KS1*J2(M+1)  > 

G(3,l)=-( (0.5*KS1**2-M*(M-1))  * Ji (M) -2 .*KL1*J1 (M*l)) 
1*I*((0.5*KS1**2-M*(M-1))*Y1(M)-2.*KL1*Y1 (M*l) ) 

G(3,2)=-M*(M*1)*(  (M-1)*J2(M)-KS1*J2(M*1)  > *I*M* (M*l) * ( (M-l ) *Y2 ( M) 
1-KS1*Y2(M*1) ) 

KS=KS1**2/KS2**2 

R=RH02/RH01 

G(3,3)*R*KS*((0.5*KSZ**2-M*(M-1>  >*J3(M>-2.*KL2*J3(M»1>) 
G(3,4)=R*KS*M*(M*1) *< (M-i)  *J4 I  Ml -KS2*  J4  <  M ♦ 1 1 ) 

G(3,5)=CIM*( (0.5*KS1**2-M*(M-1))*J1(M)-2.*KL1*J1(M*1) ) 
1*SIM*(M*(M+1)  *(  (N-1)*J2(H)-KS1*J2(M*1))) 

G(4,l)=-( IM-1I*J1 (M)-KL1*J1(M*1)>*I*  C ( M- 1 )  *  Y 1 (M) -KLl’Yl (1*1) ) 
G(4,2)*-( (1.-M**2*0.5*KS1**2)*J2(M)-KS1*J2(M*1) )  ^IM  (1-M**2 
1 ♦0«5*KS1**2)  *Y2(M)-KS1*Y2(M*1) » 

G(4,3)=R*KS*( (M-1)*J3(M) -KL2*J3(M*1) > 

G(4,4)=R*KS*( (l.-M**2+0.5*KS2**2)*J4(M)-KS2*J4(M*l>) 

G(4,5)=CIM*( (H-l)*Jl(M)-iai*Jl(M*l))*SIM*(  (l.-M**2*0.5*KSl**2) 
1*J4(N)-KS2*J4(M*1)> 

C  THE  AUXILLARY  MATRIX  FOR  THE  SOLUTION  BY  THE  CROUT  METHOO  FOLLOWS 

00  114  MXsl ,  N1 

114  A(MX,1)=G(MX,1) 

00  121  NX=2,N1 
00  121  MX=2 « N2 
A(l,MX)=G(i,MX)/A(l,l) 

SUMPRsCMFLX (  0  •  ,  0  «  ) 

IFINX-MX) 115,117, 119  ' 

115  LL*NX-1 

00  116  MM  =  1  ,LL 

116  SUHPRaSUMPRtA(NX,MM)*A (MM, MX I 
A(NX,MX)x(G(NX,MX)-SUMPR)/« (NX, NX) 

GO  TO  121 

117  LL*NX*1 

00  118  HMsl,LL 

118  SUMPR«SUMPR»A(NX,MM) *A(HM,KX> 

A (NX, NX) >G( NX, NX) -SUM PR 

GO  TO  121 

119  LL*MX- 1 

00  120  MM»1,LL 

120  SUMPR«SUMPR»A(NX,MN)*A(MM,MX) 

A (NX, MX) aG(NX,MX)-SUMPR 

121  CONTINUE 

C  THE  FOLLCNING  DETERMINES  THE  FINAL  MATRIX  BY  THE  CROUT  METHOD 
F (Nl)aA(Nl «N2) 

LL»N1-1 


c 
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00  125  HX=1»  LL 

L»N1-HX 

ML»L*1 

SUHPR*CHPLX(0.,0.) 

122  IF(Nl-ML) 124,123,123 

123  SUMPR=SUMPR»A(L,MU*F(HL) 

HL*ML*1 

GO  TO  122 

124  F(L)*A(L.N2)-SUNPR 

125  CONTINUE 

00  126  L*1,N1 

126  KK(L)*F(L)*CCNJG(F(L>> 

GAHNA*GAMMAM SIM/FI  >M2*H»1I*  C  (KS1/KL 1)  *KK  (  1>  *M*  (M*  1 )  *KK  ( 2 )  > /CS1 
1MCIM/PI)  *(2*H*1)MKK  (1)  ♦ (KL1/KS1)  *N*(M*1>  *KK(2> )/CLl 

127  CONTINUE 
GAHMA*8.686*GAMHA 
RETURN 

END 

SUBROUTINE  UPRLIF 
IMPLICIT  REAL  (J,K> 

OINENSION  NN  (4 ) 

COMMON  KL1,KS1,KL2«KS2, RH01 ,RH02, SIM, CIM.G AMMA , AK , J (4001, NMAX, 

IV (400) ,CL1,CS1,CL2,CS2 

THIS  SUBROUTINE  SETS  THE  UPPER  LIMIT  SUCH  THAT  ALL  SPHERICAL 
BESSEL  FUNCTIONS  CF  HIGHER  ORDER  CAN  BE  APPROXIMATED  TO  ZERO. 

00  105  11-1,4 

GO  TO  (96,97,98,99) ,11 

96  AK-KL1 

GO  TO  100 

97  AK=KS1 

GO  TO  100 

98  AK=KL2 

GO  TO  100 

99  AK>KS2 

100  NAK=IFIX (AK) 

N*NAK+10 

00  101  IsN, 1 950 , 5 
L*I 

FI*FLOAT (I) 

SECA*AK/(FI»0.5) 

TACA»SQRT(1.-SECA**2) 

COSA*l./SECA 

SINA*SQRT(1.*C0SA**2) 

ALP»AL0G(C0SA*SQRT(C0SA**2-l. ) ) 

0ELT*EXP ((FI*0.5)*(TACA-ALP))/(2.*(FI*0.5) *SQRT (SIN A) ) 
IFCDELT-l.E-27) 102,102,101 

101  CONTINUE 
GO  TO  107 

102  N*L 
FN«FLOAT(N) 

NP«N*1 

NT«N*2 

NU«N»10 

EXPRsEXPd.) 

00  103  I*NU, 1890  ,5 
L*I 

FI*FLOAT(I) 

0»Q. 434294 


i 

i 

I 
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Al*2.*(FI-FN*i.)*C*ALCG(AK) 

A2*(2.*FI-2.*FN*1.)*0*ALOG(EXPR) 

A3*(FN*2.5)*0*ALCGCFN*2. I ♦ (FN- 1.5 > *0* ALOG (FN-2 . > 
A4=IFI*3.5)*0*ALCG(FI*3.> 

A5*(FI-0.5) *0*ALCG  €  F I  —  1 .  I 
E* A1+A2+ A3-A4-A5 
IFC£*10.) 104,104,103 

103  CONTINUE 
GO  TO  107 

10*  NN(II)=NT 

105  CONTINUE 

NMAX=MAX0 I NN ( 1) , NN ( 2 1 ,NN ( 3 > ,NN (4) ) 

NMAX  IS  THE  UPPER  LIMIT  OF  THE  SUM  IN  THE  SCATTERING  CROSS 
SECTION.  ALL  HIGHER  TERMS  CAN  BE  CONSIOEREO  TO  BE  ZERO. 

106  RETURN 

107  MRITE(6,108) AK 

108  FORMAT  (26H  ARG  BEE  FTN  TOO  LARGE  KA= .E20.8t5X.14H  END  EXECUTION) 
STOP 

ENO 

SUBROUTINE  SPHBES 
IMPLICIT  REAL  (J.K) 

DIMENSION  R(400) *RJ (400) 

COMMON  KL1, KS1, KL2.KS2.RHOl.RHC2, SIM, CIM, GAMMA, AK,J  14  00) .NMAX. 

1Y (400) .CL1.CS1.CL2.CS2 

THIS  ROUTINE  CALCULATES  THE  SPHERICAL  BESSEL  FUNCTION  (J(I))  OF 
ARGUMENT  AK  TO  ORDER  NN.  THE  SOLUTION  IS  VALID  FOR  ARGUMENTS  AS 
LON  AS  0.05. 

NDIM=400 
00  100  1*1.400 
100  J(I)=0. 

104  NU=NMAX*10 
NT = NMAX* 2 
NP=NMAX*1 
IJ=NU 

R (NU*1) *0  . 

00  105  1=1, NU 
L*I 

SJ=FLO AT ( I J) 

R(IJ)=AK/(1.*2.* SJ- AK*fi ( I J *  1 )  ) 

IF <R (I J)-l.) 105,105,106 

105  IJ=IJ-1 
GO  TO  107 

106  IF(IJ-2) 107,107, 111 

107  RJ ( NU*1 1)  =R ( NU) 

RJ (NU) =1 . 

IJ*NU 

NUP*NU-1 

00  108  1*1, NUP 
IJ=IJ-1 
SJ*FLO AT ( I J) 

RJ(IJ)*(l.+2.*SJ)/AK*RJ(lJ*l) -RJCIJ*2) 

108  CONTINUE 

ALPM* (R J ( 1) -AK#RJ (2))*C0S(AK) ♦ AK*S IN ( AK) *R J ( 1 > 

00  109  1*1, NT 

109  J( I) *RJ  < I ) / ALPM 

110  CONTINUE 
RETURN 


o  n  o 


111  RJ(IJ+1)*R(IJ) 

LAM2=IJ*2 

IF (LAM2-NP) 112,115,119 

112  RJ(IJ) *1. 

IJsIJ-1 

l*L»l 

00  113  I*L,NU 
SJ=FL0AT  (I J) 

RJ(IJ)=(1.+2.*SJ)/AK*RJ(IJ+1) -R J ( I J+2 ) 

113  IJsIJ-1 

00  114  I=LAM2 , NT 

114  RJ(I>=RJ(I-1) »R< I- 1 > 

ALPH=(RJ(l)-AK*RJ(2n*C0S(AK) ♦AK*SIN(AK>  *RJ(1) 

N0M*NDIM-1 

IF (NT- NON 1 116,116,115 

115  NT=NOM 

116  00  117  1*1, NT 

117  J(I)*RJ( II/ALPH 

118  CONTINUE 
RETURN 

119  HRITE(6, 120) 

120  FORMAT  C25H  LAMB0A-»2.GE.N»i  END  XEQ) 

GO  TO  116 

END 

SUBROUTINE  SPHNEU 

COMMON  KL1, KS1, KL2, KS2, RH01.RH02, SIN, CIM, GAMMA, AK,J(400) ,NMAX, 

1 Y (400 ) , CL 1, CS1.CL2.CS 2 

THIS  ROUTINE  CALCULATES  THE  SPHERICAL  NEUMAN  FUNCTIGMYdl) 

OF  ARGUMENT  AK  AKC  UP  TO  CROER  NMAX.  THE  VALLES  ARE  CALCULATED 
BY  UPWARD  RECURSION  RELATIONS. 

ZERO=0. 

00  100  1-1,400 

100  Y(I)<0. 

Y (1>*ZERC-C0S (AK) /AK 

Y(2)*ZER0-C0S(AK)/AK**2-SIN(AK)/AK 

NT* NMAX ♦ 2 

00  101  1*3, NT 

AJ*FLOAT ( 1 1 -2  . 

Y(I»*( ( 2 .* AJ ♦ 1 . ) / AK ) * Y (I-l)-Y (1-2) 

101  CONTINUE 
RETURN 
ENO 
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